Course Project (MMG640 / MVE080)

Visualization of outer solar system simulation using blender

Name: Fatemeh Alibakhshi

The general process of the project is based on the Python code below.

import numpy as np
import bpy

# Gravitational constant when the mass of the sun is 1.
G = 2.95912208286e-4

# Planet names and order
planets = ('Sun','Jupiter','Saturn','Uranus','Neptune','Pluto')

# The data below is obtained from here: https://ssd.jpl.nasa.gov/horizons.cgi

# Masses relative to the sun (the increased sun mass is to account for the inner planets)
masses = np.array([1.00000597682, 
                   0.000954786104043, 
                   0.000285583733151, 
                   0.0000437273164546, 
                   0.0000517759138449, 
                   6.571141277023631e-09])

# Positions of the planets in astronomical units (au) on September 5, 1994, at 0h00 GST.
positions = np.array([[0., 0., 0.],
                    [-3.502576677887171E+00, -4.111754751605156E+00,  9.546986420486078E-02],
                    [9.075323064717326E+00, -3.443060859273154E+00, -3.008002285860299E-01],
                    [8.309900066449559E+00, -1.782348877489204E+01, -1.738826162402036E-01],
                    [1.147049510166812E+01, -2.790203169301273E+01,  3.102324955757055E-01],
                    [-1.553841709421204E+01, -2.440295115792555E+01,  7.105854443660053E+00]])

# Velocities of the planets relative to the sun in au/day on September 5, 1994, at 0h00 GST.
velocities = np.array([[0., 0., 0.],
                    [5.647185685991568E-03, -4.540768024044625E-03, -1.077097723549840E-04],
                    [1.677252496875353E-03,  5.205044578906008E-03, -1.577215019146763E-04],
                    [3.535508197097127E-03,  1.479452678720917E-03, -4.019422185567764E-05],
                    [2.882592399188369E-03,  1.211095412047072E-03, -9.118527716949448E-05],
                    [2.754640676017983E-03, -2.105690992946069E-03, -5.607958889969929E-04]])

# Compute total linear momentum
ptot = (masses[:,np.newaxis]*velocities).sum(axis=0)

# Recompute velocities relative to the center of mass
velocities -= ptot/masses.sum()

# Linear momenta of the planets: p = m*v
momenta = masses[:,np.newaxis]*velocities

# Function for Newtonian acceleration field
def acc(x, masses = masses, G = G):
    N = masses.shape[0]
    d = x.shape[-1]
    dx_pairs = x[:, np.newaxis] - x[np.newaxis, :]
    msq_pairs = masses[:, np.newaxis]*masses[np.newaxis, :]
    
    # Remove self-self interactions
    dx_pairs = np.delete(dx_pairs.reshape((N*N,d)),slice(None,None,N+1), axis = 0).reshape((N,N-1,d))
    msq_pairs = np.delete(msq_pairs.reshape((N*N)),slice(None,None,N+1), axis = 0).reshape((N,N-1))
    
    # Compute pairwise distances
    dist_pairs = np.sqrt((dx_pairs**2).sum(axis=-1))
    
    # Compute the gravitational force using Newton's law
    forces = -G*(dx_pairs*msq_pairs[:,:,np.newaxis]/dist_pairs[:,:,np.newaxis]**3).sum(axis=1)
    
    # Return accelerations
    return forces/masses[:,np.newaxis]

# Select time step and total integration time (measured in days)
h = 100 # Time stepsize in days
Frames = 365
totaltime = 100*Frames # Total simulation time in days

# Preallocate output vectors at each step
t_out = np.arange(0.,totaltime,h)
x_out = np.zeros(t_out.shape + positions.shape, dtype=float)
x_out[0,:,:] = positions
v_out = np.zeros_like(x_out)
v_out[0,:,:] = velocities

# Use Symplectic Euler method for integration
for x0, x1, v0, v1 in zip(x_out[:-1],x_out[1:],v_out[:-1],v_out[1:]):
    x1[:,:] = x0 + h*v0
    v1[:,:] = v0 + h*acc(x1)

# -------------------------
# Add the Blender code here


#set end frame of sequence
bpy.context.scene.frame_end = 365

#create planets
def Cr_UvSphere(List_Names):
    for i in List_Names:
        bpy.ops.mesh.primitive_uv_sphere_add()
        bpy.context.object.name = i
        bpy.ops.object.shade_smooth()
        
#add Scale to planets
def Scale():
    for i, item in enumerate(planets):
        myobj = bpy.data.objects[item]
        myobj.scale.x = masses[i]
        myobj.scale.y = masses[i]
        myobj.scale.z = masses[i]
        
# add movement and keyframes to planets        
def Movement():
    for i, item in enumerate(planets):
        myobj = bpy.data.objects[item]
        for f in range(Frames):
            myobj.location.x = x_out[f,i,0]
            myobj.location.y = x_out[f,i,1]
            myobj.location.z = x_out[f,i,2]
            myobj.keyframe_insert(data_path='location', frame=f+1)



Cr_UvSphere(planets)
#Scale()
Movement()

The given code stores the coordinates of the planets inside the x_out variable.

There are three functions defined here. The first function is Cr_Uvsphere, which places the planets in the origin coordinates and creates them based on the specified names. The second function is called Scale, which I have not used. Because the scale determines the size and smallness of the planets, and if it was used, the planets would become very small compared to the sun and could not be seen. The third function is Movement, which determines the movement of the planets. The total output is inside the x_out variable, which depends on three inputs: frame number, planet number, and location. After loops are defined for all planets and frames to be able to do the animation in total.

After the Python code, I downloaded the textures based on the given site and set them on the objects. To set, you first select the object and refer to the material properties and create a new material using this section. For general lighting, I downloaded an HDRI photo and went to the world properties section and set the color to the environment texture, and opened the downloaded photo in this section. Regarding the overall lighting of the scene, no light was used and only the same HDRI photos were used, which have a high dynamic range and the range of their color numbers is more than one, so they are very suitable for the lighting of the scene.

Perhaps it would be better to talk a little more specifically about sunlight. For this purpose I used a spherical point light, I put it instead of the sun but slightly bigger than it so that it can make the planets that rotate towards it brighter than behind them. Of course, it should be noted that the emission given to the material of the sun was a little luminous, but it was not able to illuminate the environment. Therefore, it had to be given light to make the space brighter.

One of the challenging parts of the project is creating orbits that are shown as a line based on the motion of the planets. For this purpose, a motion path can be used. After making the animation of the planets, the motion path can be activated. Then go to the scripting section and use the following script in new:

## path from -calculated- motion path

import bpy
ob = bpy.context.object
mp = ob.motion_path

if mp:
    path = bpy.data.curves.new('path','CURVE')
    curve = bpy.data.objects.new('Curve',path)
    bpy.context.scene.collection.objects.link(curve)
    path.dimensions = '3D'
    spline = path.splines.new('BEZIER')
    spline.bezier_points.add(len(mp.points)-1)
    
    for i,o in enumerate(spline.bezier_points):
        o.co = mp.points[i].co
        o.handle_right_type = 'AUTO'
        o.handle_left_type = 'AUTO'

After running the script, a curve is created. Returning to the layout again and entering the curve option from the settings section, the circuit can be created by setting the depth in the geometry section. In the end, the work is done, with the EV rendering engine.

I have also used this site to achieve my goal and the script that helped me most of the time. https://blenderartists.org/t/convert-object-motion-path-to-curve/679316/2

Rendered Video Using EEVEE Render Engine:

from IPython.display import Video
Video('C:\\Users\\acer\\Desktop\\Solar System.mp4', html_attributes="controls width=100%", embed=True)